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We calculate the nonlinear cotunneling conductance through interacting quantum dot systems in 
the deep Coulomb blockade regime using a rate equation approach based on the T-matrix formalism, 
which shows in the concerned regions very good agreement with a generalized master equation 
approach. Our focus is on inelastic cotunneling in systems with weakly broken degeneracies, such 
as complex quantum dots or molecules. We find for these systems a characteristic gate dependence 
of the non-equilibrium cotunneling conductance. While on one side of a Coulomb diamond the 
conductance decreases after the inelastic cotunneling threshold towards its saturation value, on the 
other side it increases monotonously even after the threshold. We show that this behavior originates 
from an asymmetric gate voltage dependence of the effective cotunneling amplitudes. 

PACS numbers: 85.35.Be, 73.63.kv 



I. INTRODUCTION 

Quantum dot devices, or so-called artificial atoms, con- 
sist of a small electronic nanostructure tunnel coupled to 
source and drain leads. In the Coulomb-blockade regime, 
sequential (one-electron) tunneling transport is exponen- 
tially suppressed and processes where two or more elec- 
trons tunnel simultaneously become the dominant trans- 
port mechanismi Among such correlated tunneling pro- 
cesses, cotunneling has received a lot of interest in re- 
cent years. Cotunneling denotes a two-electron tunneling 
process which can transfer an electron coherently from 
source to drain by a virtual population of an energet- 
ically forbidden charge state of the nanostructure. As 
energy is gained from the voltage drop during the elec- 
tron transfer, a cotunneling event can leave the structure 
in an excited state, in which case one speaks of inelastic 
cotunneling. Otherwise the energy state of the island is 
left unchanged and the process is called elastic. 

Inelastic cotunneling spectroscopy has turned out to 
be a useful tool to identify electronic, magnetic and vi- 
brational excitations in semiconducting 2,3 or carbon nan- 
otube based^ - — quantum dots as well as in single-molecule 
j unctions £ — Most importantly, the positions of conduc- 
tance peaks provide a very direct fingerprint of the ex- 
citation spectrum of the tunnel-coupled nanostructure, 
but also the more detailed bias-dependence or line-shape 
of such inelastic cotunneling peaks contains valuable in- 
formation. By now, it is well understood^— how the 
non-equilibrium pumping of excited states by the applied 
bias voltage can give rise to a cusp in the region where 
the bias voltage matches the relevant excitation energy. 
This effect is maximal for a symmetric setup. Having 
very different tunnel couplings to respectively source and 
drain electrodes implies that the nanostructure (dot or 
molecule) is almost equilibrated with one electrode and 
this effect no longer shows up. This was confirmed by an 
experiment by Parks et al. where the opening and clos- 
ing of a mechanical break junction holding a Cqq molecule 



was shown to correlate with the weakening and strength- 
ening of such non-equilibrium cusps near the threshold 
for excitation of a vibrational mode in the system^ 

In the present paper we investigate the non-equilibrium 
cotunneling in a variety of complex quantum dot sys- 
tems. For systems designed to have low-energy excita- 
tions arising from weakly broken degeneracies we find 
a characteristic gate dependence of the non-equilibrium 
modulation of the inelastic cotunneling steps. On one 
side of the Coulomb diamond, we find the characteristic 
cusped increase in conductance at threshold, but on the 
other side of the diamond this turns into a weakening of 
the conductance at threshold which renders the nonlin- 
ear conductance entirely monotonous in bias. We show 
how this comes about by an asymmetry in the cotunnel- 
ing amplitudes for processes which add or remove one 
electron from the dot ("particle- hole asymmetry"). This 
adds an important piece of information to the spectro- 
scopic toolbox insofar as such a non-equilibrium depres- 
sion of the cotunneling step in a symmetrically coupled 
device should not be mistaken for a tunnel broadened or 
thermally smeared cotunneling step in a very asymmetri- 
cally coupled system: such causes would lead to depres- 
sion of the cotunneling steps at both sides of the Coulomb 
blockade diamond. In contrast, the pecularity of our in- 
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FIG. 1: Sketch of the setups for different quantum dot systems 
investigated in this work: a double dot (DD), a triangular dot 
(TD), and a "benzene" quantum dot. 
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trinsic effect is that a depression will occur only on one 
side, while on the other the conductance will retain the 
typical cusped increase. 

The systems we study are all coupled symmetrically to 
source and drain so as to maximize the non-equilibrium 
effects under scrutiny. They are chosen according to in- 
creasing complexity, namely a lateral double-dot (DD), 
a triangular triple-dot (TD) and a benzene molecule. In 
all cases only two dots (or sites in the benzene) are tun- 
nel coupled to the leads, see Figure HJ For the triple dot 
as well as for benzene, this induces not only a breaking 
of the symmetry under rotations by n • 120° respectively 
n • 60° (n G Z), but also a degeneracy lifting: While the 
on-site energies of all uncoupled sites can be normalized 
to zero, the contact sites have to be endowed with a dif- 
ferent on-site energy £ to mimic the symmetry breaking 
which is likely to result from either tunneling renormal- 
ization 7 or electrostatic effects^ 7 -. 

To calculate the current and other observables in trans- 
port through quantum dots in the Coulomb blockade 
regime there are a number of techniques, each with their 
advantages and limitations. Following a real-time trans- 
port approach as described in RefsJ^"— , one can trace 
out the leads degrees of freedom to derive a formally ex- 
act generalized master equation (GME) for the reduced 
density matrix (RDM) of the system. This approach al- 
lows for a systematic expansion in the tunneling Hamilto- 
nian Ht, thus capturing sequential tunneling from con- 
tributions of second order in Ht and cotunneling as a 
fourth order process. The GME has the advantage that 
it is exact to a desired order. However, already at fourth 
order in Ht the number of terms is quite large and there- 
fore it is often practical to use simpler approaches which 
capture only the most relevant contributions for each or- 
der of the perturbation theory. 

In this work, we focus entirely on the Coulomb block- 
ade regime and we shall therefore calculate cotunneling 
rates within the T-matrix approach. 22 This approach is 
much simpler than the GME, and will be valid deep in- 
side a Coulomb diamond, even with the further approx- 
imations of i) neglecting the sequential tunneling contri- 
butions, and ii) approximating the denominators in the 
rates as independent of the lead electron energy. To en- 
sure that all details in the lineshapes for which we aim are 
correct, we benchmark the above mentioned approxima- 
tions as well as the T-matrix technique itself by a quanti- 
tative comparison to the GME approach. In general, we 
find the T-matrix and the GME approach to be in good 
agreement when additional effects due to level shifts and 
broadening are irrelevant, namely in the regime where 
the tunneling induced level broadening is much smaller 
than the temperature. 

The paper is organized as follows: In Section [III we in- 
troduce the model Hamiltonian and the relevant expres- 
sion to calculate the current and conductance in terms of 
transition rates is provided. In the end of this section, we 
discuss different approximations to the rates. These ap- 
proximation schemes are compared against exact fourth 



order results in Section IIIII for the case of a double-dot 
model. The triple dot and the benzene molecule are 
investigated in Section [TV] and the results are analyzed 
analytically using the simplest of these approximations. 
Conclusions are found in Section IVl 

II. THE MODEL 

For convenience the conventions e = 1 and H = 1 will 
be used throughout the paper. 

A generic quantum dot coupled to source and drain 
leads is described by the Hamiltonian 

H = Hqd + i^leads + #T- (1) 

Source and drain leads are represented by two reser- 
voirs of non-interacting electrons: i^i ea ds = Y2 a ka( €k ~ 
^a) c akcr c ak(ji wnere ol = L,R stands for the left or right 
lead. The chemical potentials fi a of the leads depend 
on the applied bias voltage which is assumed to be 
applied symmetrically across the two junctions, so that 
Ml,r = Mo =t ¥jr- In the following we will measure the 
energy starting from the equilibrium chemical potential 
fio = 0. The Hamiltonian of the quantum dot itself de- 
pends of course on the underlying nanostructure, be they 
quantum dots defined on semiconducting heterostruc- 
turesj22 carbon nanotubes 2 ^ or molecules bridging two 
contacts 25 . We consider here some archetypal model for 
quantum dots. They can be described as two or more 
(M) localized states coupled among each other. Fol- 
lowing the semi-empirical modeling of benzene j2£r— we 
study here model Hamiltonians given by 

M M 

3 = 1 cr j = l 

+ b V d iv d jv + V i3 n ™ n 3°' > ( 2 ) 

<ij> era' i<j 

where €j is the on-site energy of the level j. In Section 
IIVI we will set €j = 0, except for the two sites coupled 
to the contacts, for which we assume ej = £ <C U,Vij, 
establishing a weakly broken degeneracy. The parame- 
ter bij < describes the hopping of electrons between 
nearest neighboring states z, j. U accounts for the on- 
site charging energy, while Vij is the interaction between 
two states. The first term accounts for the influence of a 
gate voltage V g , with n being the gate coupling parame- 
ter. Its actual value strongly depends on the fabrication 
technique used for creating the quantum dot^ and typ- 
ically ranges in order of magnitude from 10 _3 -1. As it 
simply acts as a scaling factor on the gate voltage, we 
can set n = 1 throughout this paper without loss of gen- 
erality. The leads couple only to some of these localized 
states and the corresponding tunneling Hamiltonian is 
described by 

H T = J2 ( ta * d L c «k. + t a cl ka d ja<7 ) , (3) 

kaa 
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where dj aCr creates an electron in the single particle state 
\ja) which couples to lead a. The tunneling Hamiltonian 
Ht is treated as a perturbation to #qd + Pleads- In this 
work we investigate a double dot (DD), a triple dot (TD) 
and a benzene molecule, corresponding to M = 2,3,6 
in Eq. ([2]), connected to source and drain as illustrated 
schematically in Figure [TJ We give an overview of their 
relevant properties in Sections IIIII and IIVI As we shall 
show, weakly broken degeneracies in the TD and benzene 
give rise to qualitatively different inelastic cotunneling 
profiles at different gate- volt ages. This is in contrast to 
the DD which has a non-degenerate ground state and 
therefore less structure in its cotunneling amplitudes. 



A. Calculating the cotunneling current 

To make calculations easier, we shall assume T <C 
fc#T <C £ <C Ec- Here, Ec is the addition energy, 
which can in principle be expressed in terms of U and 
Vij, but with increasing number of sites in an increas- 
ingly unhandy way. The value of the level broadening 
T must stay well below the thermal energy ksT in or- 
der to justify a perturbative approach to transport. The 
degeneracy lifting £ is much smaller than the addition en- 
ergy Ec, but exceeds both T and the thermal energy by 
far. In this case, the reduced density matrix remains di- 
agonal, because one can exclude coherences between the 
no longer degenerate states, and the rate equations are 
simpler j 21 i 3Q i 31 Moreover, normal metal leads are consid- 
ered, such that a rate equation approach to transport is 
sufficient. 

For the sequential tunneling rates, a GME ap- 
proach^— yields the same result as Fermi's golden rule 
with Ht being the perturbation. The latter scheme can 
be iterated to include higher order tunneling processes 
by making use of the T-matrix 



from which transition rates from the initial to the final 
state can be calculated up to a given order in Ht- 
In general, the emerging T-matrix rates differ from 
the corresponding GME rates. The latter are exact 
to a given order in perturbation theory, and explicitly 
exclude all reducible terms j22r— i.e. divergences caused 
by the denominator in Eq. (|4]) going to zero. By 
construction the T-matrix misses in each order contri- 
butions guaranteeing these exclusion via a cancellation 
of reducible terms. Therefore unavoidable divergences 
emerge with the T-matrix technique from fourth order 
in the perturbation and onwards j 35 i 37 Meanwhile, regu- 
larization schemes to remove the divergence appearing 
in the fourth order T-matrix rates have become stan- 
dard^ 2 -"— and the T-matrix approach has been applied 
to various setups, e.g. to a double dot structure^ 3 - or 
to molecular systems where electronic and vibronic 
degrees of freedom can be strongly couplec—. In this 
context, it is important to stress that the standard 
way of regularizing the T-matrix rates does not exactly 
reproduce the GME intrinsic regularization, but the 
discrepancy between T-matrix and exact perturbation 
theory turns out to vanish deep inside the Coulomb 
blockade. The same holds for further fourth order 
contributions included by the GME which cannot 
necessarily be brought into the form of a squared matrix 
element 3 ^ and are disregarded by the T-matrix approach. 



T(E) = H T + H T 



1 



E-H + i0+ 



T. 



(4) 



We label now the states of the quantum dot with their 
particle number A 7 ", the ^-component of their spin with 77 
and an additional quantum number with I. The T-matrix 
rate for a transition between two states \N f Vr] f ) —> \Nlrj) 
of the quantum dot system is then given by^ 



\Nlr])(N'l'ri'\ 
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■Ei 



„')■ (5) 



Here the sum is over all possible initial (i) and final 
(/) states of the overall system including the leads, 
Kjvzv) = \N'l'r)')\iL)\iR), weighted by a thermal 
distribution function Wi N , lt , . The rate from Eq. (J5j) 
comprises the dominant fourth order contributions deep 
inside the Coulomb diamonds, namely the cotunneling 
effects. 

The rate equation describing the dynamics of the oc- 



cupation probabilities of the states reads 

pNi v = ^ 

- T\N'l'<n')(Nl<n\ P Nlr) + ^2 ^\Nlf])(N f l f r] f \ P N 1 V 5 

N'l'ri' N'1'r,' 

where P Nlr] (t) is the probability of finding the dot in the 
state \Nlrf) at time t. The stationary solution (t — >• 00) 
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is therefore given by 



with matrix elements 



^ F\N'l'r)')(Nlr)\P<tilt ~ ^\Nlr))(N'l'r)'\ 
N'l'r)' N'l'ri' 

with the normalization condition 

p ni v - i 

/ j ^stat — L ' 
Nlr) 



N'l'r)' 
stat i 



(7) 

(8) 



With the help of the stationary solution, we arrive at an 
approximate expression for the current up to fourth order 
in i?T? 



1 = 1 



sequential ~t~ cotunneling 7 



(9) 



h lV , = 

IL r)r)' 
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(Nlrj\d ja<T \N + ll" V ")(N + ll"rj''\4 ,jNl'r,'\ 



E 

l"r)» 



EnI't]' — EN+ll"r)" + ^k'a' + ^0+ 

{Nl V \d) aa \N - ll" V ")(N - ll" V "\d jaia ,\Nl'i) 



(14) 



Note that the effective cotunneling Hamiltonian ( [13]) 
now takes the form of a generalized Kondo, or Coqblin- 
Schrieffer model^ depending on the symmetries of the 
states | A/777). 



with the second order, 
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and fourth order contribution 
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Here, the superscripts to the rates indicate at which lead 
a the tunneling processes take place. Truncating the gen- 
eral expression (j5]) for the rate to second order, we retain 
only one tunneling event, which can either involve the 
left or the right lead. For the stationary current flow, we 
merely need to consider the balance between in- and out- 
tunneling at one of the electrodes, and we have chosen in 
Eq. (dU the left one, a = L. 

The fourth order cotunneling events transfer an elec- 
tron fully across the quantum dot, which involves two 
tunneling events at distinct leads. Therefore the resulting 
current is given by the balance between charge transfer 
from left to right (RL) and from right to left (LR). The 
cotunneling rates emerging from Eq. (|5]) can be written 



^eff 



\Nlr)){Nl'r)'\ = 27r I (f Nlr > \ H '^ l^l'r)' ) f 



xWr NlW 5(E fNlv -E lmW ), (12) 



where H^ t is given by 



#£t = E E * a * a 'Ci^x M vi4 ffCa , fcV „ (is) 



aka a'k' a' 



1. Approximation I 

To calculate the rates (fT2|) analytically, we neglect in 
a first approximation the e kot energy dependence in the 
denominators of H^ t . This is justified for small (com- 
pared to the charging energy) bias voltages, so that the 
electrons that tunnel to and from the leads have ener- 
gies around the equilibrium chemical potential and thus 
\E N ±n»r)" -E NVr)f \ > |e fca |, \e k t a/ \. Converting the sums 
over /c, k' into integrals assuming a flat band with con- 
stant density of states, a simple integration leads to the 
cotunneling rates 

^NrMNrj'l'l = 2?I " ^2 V L^r{-{E N Vt)' ~ E mrj ) - V h ) 
aa' 

X | 5g'L5gRt a t a h^^TlB (—(Ejsfl'rf — Ejsfl v ) ~ Vfc) , 
aa' 

(15) 

where ub(x) = e xp(/L)-i is the Bose function, f3 is the 
inverse temperature and v a is the density of states in lead 
a. This approximation is valid for gate and bias voltages 
inside the TV-electron Coulomb diamond. We refer to this 
approximation as Aprxl . 



2. Approximation II 

Alternatively, to get a more precise description of the 
inelastic cotunneling conductance (when hl — jiR > 
Env — Eno), we can take into account the energy de- 
pendence of H^ t . By shifting the integration variable 
£ka —> £ka + l^a in Eq. (fT4]h we see that H^ t now explic- 
itly depends on \±l^r and therefore on the bias voltage. 
In the rates, we get expressions of the form 

r ~ J def(e) (1 - f(e + fi L - fi R + E Nlr} - E Nlw )) 

(16) 
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' e - E x ± i0+ e - E 2 ± i0+ ' 



where E\ and E2 depend on lr] and Vrj' and the sum- 
mation indices in h 1 ^ If E\ = T?2, Eq. ([TB]) cannot 
be evaluated directly, because of divergences stemming 
from second order poles. This problem was stated al- 
ready in 1994^ and a regularization scheme has been 
developed and become standard within the T-matrix ap- 
proach to transport ] 32 ! 34 In this regularization scheme, a 
finite width 7 ~ T is attributed to the states which enter 
the denominators as imaginary parts. This level broad- 
ening physically stems from the tunnel coupling, but is 
not taken into account by the T-matrix approach. Thus 
the poles are shifted away from the real axis so that the 
integral can actually be performed. The resulting expres- 
sion can be expanded in powers of 7 and the leading term 
is found to be of order I/7. Together with the prefactor 
of the rates, T 2 , this term is identified to be a sequential 
tunneling term. It is excluded to avoid double counting of 
sequential tunneling processes. The next to leading order 
term is of order 7 and gives the regularized cotunneling 
rate. At this point, the actual value of the broadening 
does not matter and the limit 7—^0 can safely be taken. 
The calculation of the current with regularized cotunnel- 
ing processes and disregarding sequential tunneling rates 
(Sequential = in Eq. ©), we refer to as AprxII . 




FIG. 2: Spectrum of a double quantum dot, as described by 
Eq. ([2]) with M = 2. The gate voltage was set to V g = 20\b\, 
which corresponds to the center of the N = 2 diamond (see 
also Figures [3] and [4} . The splitting of the N = 1 states is 
Ai = —2b, the singlet-tr iplet splitting for the N = 2 states 
is A 2 = 0.5(V - U+ ^/I6b 2 + (U - V) 2 ). The interaction 
parameters are U = 20\b\, V = 10\b\. 



3. T-matrix 

Both Aprxl and AprxII are expected to fail when 
cotunneling assisted sequential tunneling processes be- 
come accessible ^ 13 ! 36 This can happen well inside the 
Coulomb diamond, when excited TV particle states are 
populated via inelastic cotunneling. Indeed, at the lines 
given by the equation =fV g ±EN±iriTENi / ri / +l^a — 0, V ^ 
(dashed lines inside the Coulomb diamond in Fig. [3]) 
the cotunneling rates become negative which leads to an 
ill-defined set of rate equations, unless we include also se- 
quential tunneling terms and allow also states with N ± 1 
to be populated. This is exactly the T-matrix approach, 
referred to as Tmat in the following. 



III. INELASTIC COTUNNELING IN A 
DOUBLE DOT 
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In this section we discuss inelastic cotunneling features 
of the simplest model described by Eq. (j2j), namely by a 
double dot system (DD). Additionally, it is used as a 
benchmark for the T-matrix approach Tmat as well as 
for Aprxl and AprxII against a calculation based on the 
GME approach. 

The spectrum of the DD system is shown in Figure [2] 
at V g = 20|6|, corresponding to the center of the N = 2 
diamond in Figures [3j HI The N = 1 states are even 
and odd combinations of electron states on the left and 
right dot with energies E e / Q = ±b. For the N = 2 states, 
we have a singlet ground state and an excited triple t 
state, separated by A 2 = Q.5(V-U+ ^I6b 2 + (U - V) 2 ). 



FIG. 3: Sketch of the stability diagram for the double quan- 
tum dot (DD) described by Eq. J2|) with M — 2 and param- 
eters U = 20|6|, V = 10\b\. Red lines indicate a transition 
between states with zero and one, green lines between states 
with one and two, blue lines between two and three electrons. 
Solid lines are for ground state to ground state transitions and 
define the Coulomb blockade regions, dashed lines involve ex- 
cited states. We have labelled the participating states by Ni. 
It indicates the ith N electron state, with associated energy 
Ei (for example, the one electron ground state is labeled with 
lo, the first excitation with li, and so on). Furthermore, the 
dotted lines indicate the onset of the inelastic cotunneling at 
Vb = Ai in the N — 1 and at Vb = A2 in the N = 2 diamond. 
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V g [|b|] 

FIG. 4: Logarithm (log 10 ) of the differential conductance I/V 6 
for the double quantum dot (DD) as a function of gate and 
bias voltage calculated with the T-matrix approach. One rec- 
ognizes the features in the dl/dV discussed schematically in 
Figure [3] Vertical dashed lines indicate the cuts used in Fig- 
ures [5] and [6] Parameters are as in Figures HE] together with 
k B T = 0.02|6|, T L = T R = 0.008|6|. Here, as well as in all 
following plots, the differential conductance is measured in 
units of $ x e 2 /h, with the scaling factor $ := Y L Y R /\b\ 2 . 
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V b [|b|] 

FIG. 5: Differential conductance I/Y& as a function of bias 
voltage calculated with the different approximation schemes 
discussed in Section [TT] as well as with the GME at the center 
of the N = 1 diamond corresponding to V g = 4.8 \b\. AprxII 
yields divergences in the conductance at resonances. The 
Tmat and GME show features at these positions but are well 
behaved. They agree almost exactly. Parameters are as in 
Figure [4] 
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FIG. 6: Differential conductance I/Yb as a function of bias 
voltage calculated with different approximation schemes in 
the N = 2 diamond at V g = 13\b\ (upper set of lines) and at 
the center of the diamond corresponding to V g — 20\b\ (lower 
set of lines). Due to the particle- hole symmetry of the DD, 
the two cuts at V g = 13|6| and V g = 27\b\ give exactly the 
same result. 

Notice the particle-hole symmetry of this system, which 
is responsible for the symmetry of the stability diagram 
around this value of the gate voltage. 

In Figure O we show a sketch of the stability diagram 
for the DD together with additional excitation lines. All 
the lines follow from energetical considerations involving 
the spectrum, see Figure [2j and the chemical potential of 
the leads. We focus on the energy range relevant to the 
case where the dot is singly or doubly occupied, i.e. to 
the Coulomb diamonds with N = 1 and N = 2. Red lines 
indicate positions for transitions between states with zero 
and one, green lines between states with one and two, and 
blue lines between states with two and three electrons, 
respectively. Solid lines are for ground state to ground 
state transitions and define the Coulomb blockade regions 
with N = 1 and N = 2, dashed lines involve excited 
states. The dotted horizontal lines indicate the onset of 
the inelastic cotunneling in the two diamonds. 

In Figure [H the conductance through the DD calcu- 
lated with Tmat is plotted on a logarithmic color scale. 
One can see nicely the features in the I/V 6 at the posi- 
tions of the lines in Figure [3] and the general resemblance 
of the two figures. Inside the diamonds, at Vb = Ai or 
Vb = A2 the threshold for inelastic cotunneling is visible 
as horizontal lines. The onset of the cotunneling assisted 
sequential tunneling (see dashed lines in Figure [3]) can 
be noticed best in the N = 1 diamond. Outside the 
diamonds, all the sequential tunneling lines can be seen. 

We are now going to compare the different approxi- 
mation schemes for the cotunneling rates as discussed in 
the previous section with the exact perturbation theory 
(GME). In Figure [5j we show the differential conduc- 
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tance as a function of the bias voltage at the center of 
the N = 1 diamond, as indicated by the dashed white line 
in Figure HJ We see that Aprxl yields good agreement 
with the GME only at small bias voltages. In particu- 
lar, the lineshape at the inelastic cotunneling threshold 
is not reproduced correctly, because the condition of va- 
lidity for Aprxl , A <C Ec, is not fulfilled here. We see 
that the other approaches predict an increasing differen- 
tial conductance for larger bias voltages, which can be 
understood from the e& dependence in the denominators 
in Eq. (CG~ 



AprxII and GME agree nicely as long as gate and bias 
voltages are such that one is in the innermost diamond 
defined by the dashed lines (see Figure [3]). Outside of 
this region, AprxII is no longer valid: Once the cotun- 
neling assisted sequential tunneling sets in, the cotun- 
neling rates in AprxII can become negative and the rate 
equations are ill-defined. 

Inside the overall Coulomb diamond, the Tmat and 
GME yield almost exactly the same result. Small relative 
deviations (few per cent) between the two approaches can 
be seen at the resonant lines (see inset in Figure [5]), which 
can be attributed to a certain class of terms in the rates 
not taken into account by the T-matrix^i 

In the TV = 2-particle diamond, a better separation 
of the energy scales defined by the addition energy and 
the inelastic cotunneling threshold is given. As expected, 
all approximation schemes and the GME give almost ex- 
actly the same result at the center of the diamond (see 
lower set of lines in Figure [6]). More towards the N = 1, 2 
charge degeneracy point, at V g = 13|b|, we see that Aprxl 
gives still a good qualitative description of the lineshape 
of the conductance, but as in the N = 1 diamonds it fails 
to reproduce the increase of the conductance due to the 
bias dependence of the denominators of the rates. The 
decrease of the conductance after the inelastic cotunnel- 
ing threshold is due to the non-equilibrium redistribu- 
tion of the population of the excited state. At low bias, 
only the ground state is populated, and only cotunnel- 
ing processes that do not change the occupation of the 
ground state are possible. When the bias is large enough 
to populate the excited state, the conductance suddenly 
increases due to the new possibilities of transferring elec- 
trons from left to right lead. The excited state starts 
to acquire a finite non-equilibrium population from this 
point on, and together with the increasing depopulation 
of the ground state this leads typically to a decrease of 
the differential conductance after the sudden increase at 
the threshold. We will discuss the behavior of the con- 
ductance at the inelastic cotunneling threshold in more 
detail in the next section. 

Since the DD is particle-hole symmetric, cuts through 
the N = 2 diamond at the same distances from the center 
towards the TV = 1,2 and N = 2,3 charge degeneracy 
points give exactly the same result. 
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FIG. 7: Spectrum of the triple dot at V g — 1.14|6| (center of 
the N = 3 Coulomb diamond). The addition energy Ec = 
Uff= 3 for TV = 3 can be read off as the distance between the 
N = 3 and the N — 2 ground state plus the distance between 
the N = 4 and the N = 3 ground state. The splitting A 
of the N = 3 ground state is about a hundred times smaller 
than the addition energy (see inset). Parameters are U = 5|6|, 
V = 2\bl£ = -0.1\b\. " 



IV. INELASTIC COTUNNELING IN DOTS 
WITH WEAKLY BROKEN DEGENERACIES 

Systems with slightly broken symmetries, e.g. 
molecules in a single- molecule junction, exhibit weakly 
broken degeneracies, where the splitting of the originally 
degenerate states is much smaller than the addition en- 
ergy. These systems provide a separation of energy scales 
which allows us to investigate inelastic cotunneling effects 
that are largely unaffected by charge fluctuations. From 
a technical point of view, this brings us into the validity 
range of the simplest approximation on the cotunneling 
rates (Aprxl ). 

We assume in the following a site independent hopping 
bij = b < for nearest neighbors z, j and a shift of the 
on-site energies of the contacted sites by £ = — 0.1 16|. 



A. Triangular triple dot 

A triple quantum dot (TD) system is described by Hqd 
in Eq. (j2j) with M = 3. We assume that the left lead is 
coupled to dot 1 and the right lead to dot 2. The cou- 
pling of dots 1 and 2 to the leads breaks the symmetry 
of the isolated molecule, and in accordance with our pre- 
vious statements we set thus 62 = e\ = £, 63 = 0. The 
corresponding energy spectrum is shown in Figure [7] as a 
function of the number N of electrons in the triple quan- 
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FIG. 8: Color coded cotunneling conductance of a TD well 
inside the N = 3 Coulomb diamond. The conductance in 
general is lowest at the center of the diamond, and increasing 
when moving the gate towards the diamonds with N± 1. The 
onset of inelastic cotunneling at Vb = A is clearly visible as a 
jump in the conductance. The charge degeneracy points are 
at V Q = -1.68|6| (with N = 2) and at V g = 3.96|6| (N = 



4). 
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Parameters as in Figure with additionally ksT = 8 



8 x 10- 5 |6|. 



turn dot. The gate voltage is chosen such that the lowest 
energy occurs when the TD is filled with N = 3 electrons. 

For £ = the N = 3 ground state is both spin 
and orbit ally degenerate. We label these orbit als with 
I = 1,2. For finite £, they will split by an energy of 
\E(l = 2) - E(l = 1)| = A(0, where A(f = -0.1|6|) is 
much smaller than the addition energy = ^1^3 ( see 
Figure [7J inset). The next excited state with N = 3 is 
separated by an energy comparable to the addition en- 
ergy and can thus be disregarded. 

In Figure [HJ we focus now on the situation when the 
system is filled with three electrons. For low tempera- 
tures, sequential tunneling is exponentially suppressed at 
small bias voltages, and the current is dominated by co- 
tunneling events. We show the cotunneling conductance 
calculated with Aprxl as a function of gate and bias volt- 
age. The inelastic cotunneling threshold is clearly seen 
as a horizontal line at V& = A. In Figure [9l we show 
three cuts of the cotunneling conductance at different 
gate voltages (calculated with Aprxl and AprxII ), one 
at the center, two towards the corners of the N = 3 di- 
amond, as indicated by the dashed white lines in Figure 
[HI As one can see from the comparison of the three cuts 
in Figure [H the magnitude of the conductance as well 
as the exact lineshape now depends strongly on the gate 
voltage. At the center of the diamond (at V g « 0.56 16|), 
and even better pronounced at lower gate voltages (e.g. 
at V g = — 0.60 1 6|), the conductance shows the expected 
behavior: it is constant below the inelastic cotunneling 



■ cut at V g = -0.60 |b| 
cut at V g = 0.56 |b| 

■ cut at V q = 1 .82 |b| 
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FIG. 9: Cotunneling conductance of a TD vs bias voltage at 
the positions indicated by the white dashed lines in Figure [8] 
The lineshape as well as the magnitude depend on the gate 
voltage. The solid/dashed/dot-dashed curves are calculated 
with Aprxl , their dotted companions with AprxII . The two 
approaches agree well for low bias voltages Vb <C Uff= 3 . 



threshold, and above it shows a step with a cusp. The ori- 
gin of this cusp lies in the non-equilibrium redistribution 
of the occupation probabilities of the two orbitals. For 
bias voltages below the threshold, only the ground state 
is populated. Above, the occupation probability of the 
excited state rises and increases with the bias, heading 
towards its saturation value. Until the saturation value 
is reached, the conductance will change with the bias. 
This also true for the cut at V g « 1.82 however, there 
is no cusp, but a steady further increase in conductance 
above the inelastic cotunneling step. To understand this 
different behavior, we analyze the expression for the co- 
tunneling current and the underlying rates. 

We allow only the orbitals of the split ground state 
to be populated. We therefore have to solve the rate 
equations ©, (|5J) for p N = 3lr i^ p or zero magnetic field, 
we expect p N = 31 ^ = pN=3il anc j conveniently we can 
reduce the problem to two independent variables Pi = 
T, v pN=3lri i ^ = 1,2, where the index N = 3 has been 
dropped. The following analysis is now performed under 
the assumptions that T = 0, Vb > and E(l = 1) < 
E(l = 2) without loss of generality. 

The current is then given by 



*52 T \n(i\ p i> 



IV 

and the differential conductance follows as 



dl 



1 I// 



dPi 



(17) 



(18) 



Here, rj^/w/i is the cotunneling rate for changing the quan- 
tum dot from the state I to V and thereby transferring 
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an electron from the left to the right lead. Within Aprxl 
, we can write the total rate as in ([T5]) 

l \m\ - z^Vwi 

CXOi' 

= ^2 rfn(i\ (Va'-Va + E v - Ei) (ncx'-fioc +E V -Ei), 

CXOi' 

(19) 

where 7^^| depends on the gate voltage only and the re- 
maining terms depend only on the bias voltage. We are 
especially interested in the conductance slightly above 
the inelastic cotunneling threshold, where V& = A+e, e — » 
+ . At this point, still Pi w 1 and P 2 « 0, while 



dPi 



Vb=A+e 



<0 ^ 



V b =A+e 



> and ArffZm = 



dv b A \i>)(i\ 



With these inputs, we get for the conductance 



dl 



V 5 =A+e 



1 |1>(1| 



L \m\ 



1 II 



RL 
|1)<2| 



r RL 
1 |2)<2| 



dPi 

dH 



(7|i><i| +7| 2 >< 



i|) P i + (^X2|+7|2X2|)^2. (20) 



ing bias, can win, so that the conductance does not 
show a cusp, but increases monotonously after the step. 
In other words, with a sufficiently large value of 7^ 2 |' 
the conductance will keep growing once the state / = 2 
starts to be increasingly occupied at Vb ^ A. Below 
the threshold, the elastic cotunneling conductance is set 
by the prefactor to Pi in Eq. ([20]) . At very large bias, 
however, the two states become equally populated, the 
first line in Eq. ([20]) vanishes and the with a large value 
for 7|f^ 2 |' ^ ne saturation conductance can become much 
larger than the elastic sub-threshold conductance. This 
implies that a monotonous increase of the conductance 
across the threshold will be accompanied by a large step- 
height, i.e. a large difference between the conductance 
at Vb = and at Vb ^> A. This is clearly seen to be the 
case in Figure [9] 

One can make the above statements about large cou- 
plings more precise, by inserting the stationary solutions 
for Pi and P2 given by 



Pi 



|1><2| 



|2><1| 



|1><2| 



|2><1| 



|1><2| 



|2><1| 



(21) 



From this expression, one sees that if 7^ 2 | * s l ar g e ? 
the contribution containing P2, which grows with rais- 



into the second derivative of the current at Vb 
From this, one obtains 



e. 



d 2 I 
dV* 



RL ( I RL _ RL \ I RL 

7| 2 )(i| \y(\m\ 7 i2)(iiJ \Jii)<i 



ryRL 



• 4 7|2Ki|7|^X2|) 



RL 

><2| 



. rsjLL _|_ ryRR 

7|1><2| + ~|1><2| 



(22) 



7|2Xi| (7|2Xi| + 7|?xi| " 7|?X2| " 7,2X2| i){2\ +7|i) ( 
a(2 7^2|+7^ 2| +7^ 2 |) 2 



which is positive, thus giving rise to a monotonously in- creasing cotunneling conductance, whenever 
I 



7|2)(2| 



> 



7|i)(i| 



(>y RL 

V 7 |2)(l| 



7^y(2|) (7|iX2| 



•^)(2|) 



U RL 



RL 



7| 2 xi|7|i>(2| 



syRL 



syRL 

7 |2><1| 



■ryLL 

7 |1)<2| 



syRR 

7 |1><2| 



(23) 



The question now remains, under which circumstances we have to analyze the transition amplitudes 
this condition can be fulfilled. To answer this question, 
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FIG. 10: Left panel: The quasi degenerate N = 3 ground states of a TD with / = 1, 2 as a function of the degeneracy lifting 
parameter £. For £ < 0, the state with / = 1 is the ground state, for £ > the / = 2 state has lower energy. Right panel: 
Overlap matrix elements of the levels I = 1,2 to the ground states of N =b 1 particles, a: (N — lg\ E CT dj a o-\Nl — 1), b: 
(N - lg\ d ja „\Nl = 2), c: (Nl = 1| £ ct ^ CT |iV + lp>, d: (Nl = 2| £ ct d ja<r \N + lp>. 



E 

l"7]" 



(Nlr]\d jRa \N + lZVX^V + "Vl^cr' I^V) 



Ejsii' — En 



E 



(iV^|4jiV - 1ZVX7V - \l"rjf'\d jLal \NVrjf) 
En-ii" — Eni> 



■ (24) 



We see that they depend on the overlap matrix elements 
of the tunneling Hamiltonian in the numerator, and on 
the energy differences of the states involved in the cotun- 
neling process in the denominator. A is small compared 
to the addition energy, and we keep a distance to the 
edges of the diamonds, so that for our analysis, we can 
set E N = 3 (l = 1) = E N = 3 (l = 2) in the denominator of 
the above expression. However, as we approach one of 
the two charge degeneracy points (either N o N + 1 or 
TV — 1 o TV) on the axis of the gate voltage, the contri- 
butions 



(Nl\d jRa \N + lg){N + lg\^ L JNV) 



E 



Nl' 



E 



N+lg 



or 



(Nl\d) L<T \N-lg)(N-lg\d jR JNl') 



En- 



Env 



(25) 



(26) 



with \N zb lg) being the ground states with N ± 1 elec- 
trons, are dominant in Eq. (j5]). The excited states con- 
tribute as well, but less due to the energy difference in 
the denominator, and their influence will not change the 



qualitative behavior of the differential conductance. It is 
therefore necessary to analyze separately the matrix ele- 
ments (Nl\ £ CT d ja<T \N+lg) and (N - lg\ E CT d jair \Nl'). 
These are compared in Figure [lOj where we also inves- 
tigate the behavior of the quasi degenerate levels as a 
function of the degeneracy lifting £. Plotting the ener- 
gies of the two levels with I = 1, 2 versus £, we see that 
for £ < 0, the state with I = 1 is the ground state, but 
for £ > the state with I = 2 has lower energy. The 
overall dependence of the matrix elements on £ is rather 
weak, but the coupling to the ground states with N ±1 
electrons of these two states is seen to be very different. 
The matrix element of the N — 1 ground state with I = 1 
is about twice as large as the one with I = 2, while for the 
elements with the 7V + 1 ground state, this situation is re- 
versed. This is reason why by changing the gate-voltage 
one can tune the system into a configuration where 7^ 2 | 
by far exceeds 7j^i| such that the conductance increases 
monotonously even after the inelastic cotunneling thresh- 
old. 

The cuts in Figure [9] were done for £ < 0. Choosing 
instead £ > 0, the picture would be reversed and the 
conductance at lowest gate voltages (close to the side of 
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FIG. 11: Spectrum of the benzene molecule for TV = 6 (neu- 
tral molecule), N = 7 and TV = 8 electrons. The initially 
degenerate N = 7 ground state is split due to the coupling 
to the leads (see inset). Parameters are U = 4|6|, V = 2.4|6|, 

Z = -o.i\b\. 




FIG. 12: Cotunneling conductance vs bias voltage for a ben- 
zene molecule with TV = 7 electrons coupled to the leads in 
met a configuration. The three curves correspond to three 
values of the gate voltage being closer to the TV = 6 (neu- 
tral) particle diamond, at the center of the TV — 7 diamond 
and closer to the TV = 8 diamond. Parameters are U = 4|6|, 
V — 2A\b\ (we consider nearest neighbor interaction only), 
£ = — 0.1|6|, and the remaining as in Figure [8] 



the TV — 1-diamond) would be monotonously increasing, 
while the conductance close to the TV + 1-diamond would 
now show the cusped lineshape. 



B. Benzene 

We find exactly the same effect for a singly charged 
(TV = 7) benzene molecule coupled to the leads in meta 
configuration. The spectrum of benzene exhibits a lot of 
degeneracies due to the D$h symmetry of the molecule^ 
The environment of a molecular junction can break the 
perfect symmetry of the molecule in various ways^ As 
for the triple-dot, we model this by ascribing a different 
on-site energy to the contact sites. We diagonalize #qd 
exactly, and use the eigenstates and eigenvalues for each 
charge-state to calculate all relevant cotunneling rates. 
The energy-spectrum is shown in Figure [11] and we now 
focus on the inelastic cotunneling corresponding to the 
weakly broken degeneracy in the TV = 7 state. In Fig- 
ure [121 we show three cuts through the N = 7 Coulomb 
diamond of benzene at different gate voltages, one cor- 
responding to the center and two towards the charge de- 
generacy points with TV = 6 and TV = 8. Also here, 
the lineshape at the inelastic cotunneling threshold has 
a marked dependence on the gate-voltage arising from 
a pronounced gate-voltage asymmetry in the cotunnel- 
ing amplitudes. Closer to the TV = 6 diamond, virtual 
tunneling-out processes are closer to resonance (have a 
smaller energy denominator) and closer to the TV = 8 
diamond virtual tunneling-in processes dominate. As for 
the TD, the monotonously increasing conductance closer 
to the TV = 8 diamond is clearly seen to also have a larger 
step- height. 



V. CONCLUSIONS 

In conclusion, we investigated cotunneling phenomena 
in complex quantum dot systems and demonstrated that 
systems with weakly broken degeneracies can exhibit a 
marked gate- voltage dependence of the nonlinear cotun- 
neling conductance traces. The effect relies on the non- 
equilibrium population of the excited state and is there- 
fore most pronounced in devices coupled symmetrically 
to source and drain electrodes. The inelastic cotunneling 
threshold was found to be modulated so as to become 
either cusped or monotonously increasing, depending on 
whether the strongest transport channel is via the ground 
state or via the first excited state. 

In Ref. 7 , the inelastic cotunneling thresholds were 
shown to acquire a gate-dependence due to the differ- 
ence in tunneling-induced level-shifts for the two differ- 
ent levels involved. Whereas that effect shows up with 
only tunnel coupling to a single lead, it is important to 
recognize that the gate-dependent modulation of the step 
which we discuss here relies entirely on the coupling to 
two different leads. Entering the Kondo regime for which 
such level-shifts become important, both effects could be 
observed simultaneously and therefore it would be inter- 
esting to study this stronger coupled regime more closely 
in future studies. 

From the technical point of view, we have demon- 
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strated that the widely used simplification of the T- 
matrix approach (Aprxl ) is indeed in quantitative agree- 
ment with the exact fourth order perturbation theory 
(GME) in regions of gate and bias- volt age for which se- 
quential tunneling resonances are strongly suppressed. 
With a poorer separation of energy scales, i.e. when 
the inelastic cotunneling threshold is no longer much 
smaller than the charging energy, Aprxl is insufficient 
but approximation AprxII and the T-matrix approach 
still perform very well and have a fairly large range of 
validity, within which they yield good agreement with 
the GME results. In particular, they describe very well 
the lineshape of the inelastic cotunneling conductance 
for systems with weakly broken degeneracies, i.e., with 
A <C Ec- Aprxl gives rise to substantial simplifications 



and allows writing occupation numbers and current in 
closed analytic form, despite the potential complexity 
of the underlying quantum dot systems which is now 
wrapped up in the virtual transition-amplitudes compris- 
ing the effective exchange-cotunneling matrix-elements, 
see Eq. (jHJ). As such, Aprxl can also be used for the in- 
vestigation of Kondo-enhanced inelastic cotunneling^ in 
more complex quantum dot or single- molecule systems, 
where the most relevant terms in higher order perturba- 
tion theory will be the log-singular terms underlying the 
Kondo effect. 

We thank Andrea Donarini for fruitful discussions. We 
acknowledge financial support by the DFG within the 
research programs SPP 1243 and SFB 689. 
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